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Abstract 

To assess how ecological and morphological disparity is interrelated in the 
adaptive radiation of Antarctic notothenioid fish we used patterns of opercle 
bone evolution as a model to quantify shape disparity, phylogenetic patterns of 
shape evolution, and ecological correlates in the form of stable isotope values. 
Using a sample of 25 species including representatives from four major noto- 
thenioid clades, we show that opercle shape disparity is higher in the modern 
fauna than would be expected under the neutral evolution Brownian motion 
model. Phylogenetic comparative methods indicate that opercle shape data best 
fit a model of directional selection (Ornstein-Uhlenbeck) and are least sup- 
ported by the "early burst" model of adaptive radiation. The main evolutionary 
axis of opercle shape change reflects movement from a broad and more sym- 
metrically tapered opercle to one that narrows along the distal margin, but with 
only slight shape change on the proximal margin. We find a trend in opercle 
shape change along the benthic-pelagic axis, underlining the importance of this 
axis for diversification in the notothenioid radiation. A major impetus for the 
study of adaptive radiations is to uncover generalized patterns among different 
groups, and the evolutionary patterns in opercle shape among notothenioids 
are similar to those found among other adaptive radiations (three-spined stick- 
lebacks) promoting the utility of this approach for assessing ecomorphological 
interactions on a broad scale. 
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Introduction 

Morphological disparity, a measure of the variability in 
morphological form, is well recognized to be unequally dis- 
tributed across vertebrate phylogeny (e.g., Erwin 2007; Pig- 
liucci 2008; Sidlauskas 2008). Evolutionary constraints 
place viability limits on morphological form, leaving gaps 
in phenotypic space; for instance, developmental programs 
begin at selected start points, making the achievement of 
some forms not possible along a particular ontogenetic 
pathway (e.g., Arthur 2004; Salazar-Ciudad 2006; Raff 



2007; Klingenberg 2010), and the interactions between 
genetic or phenotypic traits can channel variation in fixed 
directions (e.g., Marroig and Cheverud 2005, 2010; Brake- 
field 2006). Understanding why phenotypic spaces possess 
these properties, and the evolutionary processes underlying 
their patterning, has long captured the attention of evolu- 
tionary biologists (e.g., Wright 1932; Simpson 1953; Gould 
1989; Carroll 2005). In this regard, the study of adaptive 
radiations, groups that have rapidly diversified from a 
common ancestor to occupy a wide variety of ecological 
niches, has been of particular interest because these bursts 
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of speciation have been causally implicated in generating 
significant portions of biodiversity, or, in other words, fill- 
ing phenotypic space (e.g., Schluter, 2000; Seehausen 2007). 

Classical model examples of adaptive radiation include 
the Anolis lizards of the Caribbean (e.g., Losos 2009), 
cichlid fishes of East Africa's great lakes (e.g., Kocher 
2004; Seehausen 2006; Salzburger 2009; Santos and Salz- 
burger 2012), and Darwin's finches from the Galapagos 
(e.g., Grant and Grant 2006). These systems have been 
well studied, and thanks to a host of empirical and theo- 
retical approaches, some commonalities about the process 
of adaptive radiation have been found. All modern defini- 
tions of adaptive radiation feature a multiplication of spe- 
cies and adaptive diversification (Schluter, 2000; Gavrilets 
and Losos 2009; Glor 2010; Harmon et al. 2010). At the 
same time, however, the myriad and often lineage-specific 
interactions that guide evolutionary processes make diffi- 
cult our understanding of how well these generalities may 
fit other, less intensively studied adaptive radiations, and 
much disagreement persists regarding the meaning of 
adaptive radiation (Harder 2001; Olson and Arroyo- 
Santos 2009). A main feature of adaptive radiation mod- 
els is the idea that rapid diversification is possible under 
conditions of ecological opportunity (Schluter, 2000), and 
mathematical models predict that speciation rates and 
major ecological differences are highest at early stages of 
radiation ("early burst"), but decline as more and more 
niches become filled over time and ecological opportunity 
reduces (Gavrilets and Losos 2009). No two environments 
are the same, and the extent to which ecological condi- 
tions may place different demands on the generation and 
structuring of variation, and therefore impact our under- 
standing of adaptive radiation models, is not well known 
(Day et al. 2013). To fill these gaps, both a wider sam- 
pling of the tempo and mode of adaptive radiations and 
a focus on probing the diverse boundaries of environ- 
ments in which radiation has occurred are necessary. 

In this study we focus on the Antarctic notothenioids, a 
suborder of marine perciform fishes that represent an 
example of adaptive radiation in an extreme environmen- 
tal setting (Eastman and McCune 2000; Matschiner et al. 
2011; Rutschmann et al. 2011; Lau et al. 2012). Antarctic 
notothenioids are endemic to the Southern Ocean, the 
world's coldest and iciest marine waters (Dayton et al. 
1969; Hunt et al. 2003; Cheng et al. 2006). Together with 
the purely Antarctic Nototheniidae, Harpagiferidae, Bathy- 
draconidae, Artedidraconidae, and Channichthyidae, the 
clade also includes the three ancestral families Bovichtidae, 
Pseudaphritidae, and Eleginopidae, represented by 11 
mainly non-Antarctic species. The main radiation of the 
Antarctic group arose around 23 million years ago, near 
the Oligocene-Miocene boundary (Matschiner et al. 
2011), coincident with the development of Antarctic sea 



ice and the progressive isolation of the Antarctic shelf. In 
response to changes in water temperature, Antarctic noto- 
thenioids developed adaptive features such as antifreeze 
glycoproteins (AFGPs) and, in one family, loss of hemo- 
globin that enabled them to survive and diversify in freez- 
ing waters not habitable by other teleosts (Eastman 1993; 
Chen et al. 1997; Hofmann et al. 2005; Near et al. 2012). 
Besides their taxonomic diversity, comprising 132 pres- 
ently recognized species (Eakin et al. 2009), notothenioids 
occupy a large number of very different ecological roles 
(Eastman 1993). Several lineages independently evolved 
toward a pelagic lifestyle, a transition which, because not- 
othenioids do not possess a swim bladder, required exten- 
sive morphological and physiological adaptations to 
achieve neutral buoyancy (Klingenberg and Ekau 1996; 
Eastman 2005). The purely Antarctic notothenioids 
include five major groups that differ both in their species 
richness and extent of morphological and ecological 
diversification (Eastman 2005), these are as follows: Arte- 
didraconidae, Bathydraconidae, Channichthyidae, Harpag- 
iferidae, and Nototheniidae. The family Nototheniidae has 
undergone the most ecological and morphological diversi- 
fication, and includes 33 Antarctic species with life styles 
that range from purely benthic, epibenthic, semipelagic, 
and cryopelagic to fully pelagic (Klingenberg and Ekau 
1996; Eastman 2005). In contrast, Harpagiferidae repre- 
sents a monogeneric family of nine ecologically very simi- 
lar species, and also Artedidraconidae solely comprise 
benthic species that mainly differ in body size (Eakin et al. 
2009). Bathydraconidae are morphologically rather diverse 
and range from moderately robust to more elongate and 
delicate species, including the deepest-living notothenioids 
(DeWitt 1985) as well as shallow-living forms. Channich- 
thyids are fusiform pike-like fishes, and uniquely among 
vertebrates they lack hemoglobin. Typically living at 
depths of less than 800 m, channichthyids are quite large 
fishes (ca. 50 cm length) and most adopt a combined 
pelagic-benthic lifestyle (Eastman 2005; Kock 2005). 

Despite recent attention to the key features of the noto- 
thenioid radiation (e.g., Eastman 2005), very few studies 
have explicitly considered the evolution of morphological 
and environmental features among notothenioids (Ekau 
1991; Klingenberg and Ekau 1996), although there exist a 
large number of studies of ecomorphology and functional 
ecology for other fishes (e.g., Lauder 1983; Bemis and Lau- 
der 1986; Wainwright 1996; Westneat et al. 2005; Westneat 
2006; Grubich et al. 2008; Mehta and Wainwright 2008; 
Cooper and Westneat 2009; Holzman et al. 2012). Here, we 
collect geometric morphometric data to describe shape 
evolution for a craniofacial bone, the opercle, which articu- 
lates with the preopercle and supports the gill cover in bony 
fish. Use of geometric morphometries to analyze shape 
explicitly improves upon previous schemes of simple linear 
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measurements (Klingenberg and Ekau 1996), which may 
incur complications due to size-related effects in organisms 
such as fishes, which are characterized by indeterminate 
growth. Opercle shape is indirectly related to foraging ecol- 
ogy because besides protecting the gill cover, the opercle plays 
a primary role in the suction pump phase of the respiration 
cycle (Hughes, 1960: Anker 1974; Lauder 1979). In a simple 
distinction, fish feeding on benthic prey typically use a suc- 
tion-feeding mechanism, whereas those feeding on plank- 
tonic prey rely on ram feeding (Gerking 1994; Willacker 
et al. 2010). The ability to produce strong negative pressure 
gradients within the oral cavity is recognized as an important 
evolutionary axis of diversification (Collar and Wainwright 
2006; Westneat 2006), and additional factors such as skull 
kinesis and jaw protrusion interact in a complex way to allow 
capture of aquatic prey (Holzman and Wainwright 2009). It 
is likely that differences in opercle size and shape along the 
trophic axis affect the functionality of the suction pump. 

Using the opercle as an example of a functionally impor- 
tant and taxonomically variable craniofacial element, the 
aim of this study was to assess the interaction between ecol- 
ogy, inferred from stable isotope data, and morphology 
across the notothenioid clade, and to quantify the tempo 
and mode of ecomorphological interactions using disparity 
through time (DTT) and phylogenetic comparative meth- 
ods. Taking advantage of its relatively well-documented 
development and growth (e.g., Cubbage and Mabee 1996; 
Kimmel et al. 2005, 2008), several studies have previously 
focused on the opercle, using three-spined sticklebacks as a 
"model" system to investigate the interplay between evolu- 
tion and development. The three-spined stickleback is an 
example of a genealogically very recent species complex, 
repeatedly derived from marine ancestors after the retreat 
of the Pleistocene ice sheets to colonize freshwaters (Colosi- 
mo et al. 2005; Makinen and Merila 2008; Jones et al. 
2012a,b). Accompanying these colonizations, opercle shape 
has been shown to have repeatedly evolved along the same 
shape trajectory in geographically distinct populations, on 
a relatively short time scale, following divergence from an 
oceanic ancestor (Kimmel et al. 2008, 2011; Arif et al. 
2009). Variability in opercle shape among freshwater popu- 
lations was also found to be associated with habitat, 
differing along the benthic-limnetic axis (Arif et al. 2009). 
These results demonstrate the utility of geometric morpho- 
metries to quantify opercle shape, and imply that the glob- 
ally recovered dilation-diminution trajectory of opercle 
shape change is most likely naturally selected. Fossils are 
recognized as an important component to the study of 
adaptive radiation (Gavrilets and Losos 2009), and the op- 
ercle model further provides an opportunity to gain insight 
into the temporal persistence of evolutionary patterns of 
shape change and their implications for the paleobiology of 
extinct species flocks (Wilson et al. 2013b). 
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Material and Methods 
Sample and collection 

All specimens photographed for this study were collected 
during RV Polarstern expedition ANT-XXVIII/4 to the 
Scotia Sea in 2012. Species identification followed Gon 
and Heemstra (1990) and the FAO species identification 
sheets for fishery purposes (Fischer and Hureau 1985). 
The location, date, time, water depth, and station were 
recorded for each trawl from which fishes were photo- 
graphed (Table SI). 

The study is based on measurements of 89 specimens 
from 25 notothenioid species (Table 1, Fig. 1), including 
representatives from each of the families Nototheniidae, 
Artedidraconidae, Bathydraconidae, and Channichthyidae. 
Each specimen was photographed in a standardized man- 
ner after being fixed in position on a flat surface using large 
steel needles. A Nikon D5000 camera (Nikon Corporation, 
Tokyo, Japan) mounted on a tripod, with the camera lens 
positioned such that it was parallel to the plane of the oper- 
cle, was used to capture a close-up image of the left side of 
the head in lateral orientation. At the initial data collection 
(photography) stage, each species was represented by 



Table 1. Specimens analyzed in this study. 



Group 


Species 


N 


Lifestyle 


Bathydraconidae 


Akarotaxis nudiceps 


1 


benthic 


Bathydraconidae 


Parachaenichthys charcoti 


1 


benthic 


Artedidraconidae 


Artedidraco skottsbergi 


1 


benthic 


Artedidraconidae 


Pogonophryne scotti 


1 


benthic 


Channichthyidae 


Chaenocephalus aceratus 


3 


benthic 


Channichthyidae 


Champsocephalus gunnari 


7 


pelagic 


Channichthyidae 


Chionodraco rastrospinosus 


7 


benthic/ 








benthopelagic 


Channichthyidae 


Cryodraco antarcticus 


7 


pelagic/benthic 


Channichthyidae 


Neopagetopsis ionah 


1 


pelagic 


Channichthyidae 


Pseudochaenichthys 


3 


pelagic/ 




georgianus 




semipelagic 


Channichthyidae 


Chaenodraco wilsoni 


4 


pelagic 


Nototheniidae 


Dissostichus mawsoni 


12 


pelagic 


Nototheniidae 


Gobionotothen gibberifrons 


10 


benthic 


Nototheniidae 


Lepidonotothen larseni 


1 


semipelagic 


Nototheniidae 


Lepidonotothen nudifrons 


2 


benthic 


Nototheniidae 


Lepidonotothen 


7 


benthic 




squamifrons 






Nototheniidae 


Notothenia coriiceps 


2 


benthic 


Nototheniidae 


Notothenia rossii 


9 


semipelagic 


Nototheniidae 


Pleuragramma antarcticum 


2 


pelagic 


Nototheniidae 


Trematomus eulepidotus 


1 


epibenthic 


Nototheniidae 


Trematomus hansoni 


2 


benthic 


Nototheniidae 


Trematomus newnesi 


2 


cryopelagic 


Nototheniidae 


Trematomus scotti 


1 


benthic 


Nototheniidae 


Trematomus tokarevi 


1 


benthic 


Nototheniidae 


Trematomus bernacchii 


1 


benthic 
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LIFESTYLE 


A Channichthyidae C Artedidraconidae 
B Bathydraconidae D Nototheniidae 


0 Pelagic O Semi-pelagic O Epibenthic 
#Benthic Ocryopelagic Q Benthic/pelagic 



Figure 1. Phylogenetic relationships for the species used in this study. Filled and open circles indicate lifestyle, and major clades are highlighted 
and labeled. Phylogenetic relationships were based on those reported by Rutschmann et al. (2011) and Matschiner et al. (2011). Photographs of 
species used in this study (not to scale), from top to bottom, are as follows: Cryodraco antarcticus, Chionodraco rastrospinosus, Champsocephalus 
gunnari, Parachaenichthys charcoti, Artedidraco skottsbergi, Notothenia coriiceps, Pleuragramma antarcticum, Trematomus eulepidotus, 
Lepidonotothen squamifroms, and Dissostichus mawsoni. See Table 1 for further details of the study sample. 



between two and 30 individuals, as was available on the 
trawl, and subsequent pruning of the data set for geometric 
morphometric data collection was conducted to include 
only undamaged adult specimens, and exclude clear out- 
liers in terms of body length to minimize intraspecific allo- 
metric variation. 

Morphometric analyses 

We used an outline-based geometric morphometric 
approach to compare opercle shape across the nototheni- 
oid species examined. Geometric morphometries is a useful 
method to analyze morphological shape, capturing data 
that are easily visualized in morphospace ordinations and 
tractable to multivariate statistical methods (e.g., Book- 
stein, 1991; Adams et al. 2004; Mitteroecker and Gunz, 
2009). Here, and similar to a previous study (Wilson et al. 
2013b), an outline-based approach was chosen to assess 
interspecific shape variation because the curved nature of 
the operculum makes difficult the identification of a suffi- 
cient number of biologically meaningful, homologous, 
landmark points required for an accurate description of its 
shape across species. Eigenshape (ES) analysis is based on 



the definition of additional points of reference, or so-called 
semilandmarks (MacLeod, 1999) that are used to fill land- 
mark-depleted regions, and in doing so enable the shape 
difference located in-between landmarks to be sampled, 
and the global aspect of a boundary outline to be evaluated 
(Wilson et al. 2011). ES analysis has proven to be success- 
ful in elucidating subtle shape variation in a wide variety of 
contexts (e.g., Polly, 2003; Krieger et al. 2007; Wilson et al. 
2008; Astrop, 2011; Wilson 2013a) and is particularly suit- 
able for this study as it affords the possibility to examine 
localized variation in opercular shape. 

For each specimen, the outline of the opercle was 
traced using the software tpsDig (v. 2.16, Rohlf, 2010) 
(Fig. 2). A type II (Bookstein, 1991) landmark was 
defined as the starting point for each outline, and is 
described as the maxima of curvature on the dorsal mar- 
gin of the bone (Fig. 2). Each outline was resampled to 
create 100 equidistant landmark points. Cartesian x—y 
coordinates of these landmark points were converted into 
the phi <P form of the Zahn and Roskies (1972) shape 
function, required for ES analysis (MacLeod, 1999). ES 
analysis was performed using FORTRAN routines written 
by Norman MacLeod (NHM London). The method is 
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Figure 2. Outline-based geometric morphometric methods were 
used to capture the entire outline of the bone using 100 equidistant 
landmarks (open circles). A spatially homologous point (large color 
filled circle) was defined as starting point for each specimen. 

based on a singular value decomposition of pairwise 
covariances calculated between individual shape functions, 
and produces a series of mutually orthogonal latent shape 
vectors which represent successive smaller proportions of 
overall shape variation such that the greatest amount of 
shape variation is represented on the fewest independent 
shape axes. Each specimen has a series of eigenscores, 
representing its location along each axis, and therefore 
specimens can be projected into a multidimensional 
morphospace to visualize shape differences. Interspecific 
differences in shape were assessed using analysis of vari- 
ance (ANOVA) coupled with post hoc tests. 

Stable isotope data 

Stable isotopes of carbon and nitrogen can be used to 
provide insights into community trophic ecology because 
they show a stepwise enrichment with trophic level in 



marine systems (Hobson et al. 1994). The heavier isotope 
of nitrogen ( 5 N) is enriched by 3-4 per mil per trophic 
level and can therefore be used to infer trophic position, 
whereas the heavier isotope of carbon ( 13 C) is typically 
used to estimate the source of carbon for an organism, 
and practically applied to distinguish between near-shore 
(littoral) and open water (pelagic) environments (Post 
2002). Isotope data are expressed in delta (d) notation of 
per mil (% 0 ) versus atmospheric N 2 (AIR) and carbonate 
standards (V-PDB), using the equation S = [(R samp i e / 
Rstandard) - 1] x 1000, where R represents the ratio of the 
heavy to the light isotope (i.e., 13 C/ 12 C and 15 N/ 14 N) 
(Rutschmann et al. 2011; p4712). For all species exam- 
ined, except Akarotaxis nudiceps, Artedidraco skottsbergi, 
Trematomus scotti, and Trematomus bernacchii for which 
data were not available, stable isotope data (S 13 C and 
<5 15 N isotope) were compiled from Rutschmann et al. 
(2011) to assess the relation between opercle shape and 
lifestyle patterns. Rutschmann et al. (2011: File SI) sam- 
pled multiple specimens per species and we therefore 
computed, for each species analyzed here, an average 
value for <5 13 C and for <5 15 N. 

The relation between shape and ecology was assessed 
using phylogenetic generalized least squares (PGLS) regres- 
sion of <5 13 C with scores for axes ES1-ES8, and separately of 
<5 15 N with scores for axes ES1-ES8. PGLS uses a regression 
approach to account for phylogenetic relationships and 
assumes that residual traits are undergoing Brownian 
motion (BM) evolution (Rohlf 2001; Butler and King 2004; 
Blomberg et al. 2012). Regressions were conducted in the 
freely available statistical environment of R (http:// 
r-project.org/) using the packages "geiger" and "nlme" (gls 
function) on a pruned data set (N = 21) comprising all 
species for which we had stable isotope values. 

Disparity analyses 

To visualize the relationship between phylogeny and 
taxon spacing in ES space, phylomorphospaces were con- 
structed using ES scores. For species represented by more 
than one specimen, average scores along each axis were 
used for each phylomorphospace ordination. Following 
Sidlauskas (2008), the plot tree 2D algorithm in the rhet- 
enor module (Dyreson and Maddison 2003) of mesquite 
(Maddison and Maddison 2011) was used to construct 
phylomorphospaces for ESI versus ES2 and ESI versus 
ES3, comprising 75.4% of sample shape variance: subse- 
quent axes were not plotted as each contained less than 
8.6% of sample variance, and were not deemed significant 
under the broken-stick model (Jackson 1993). The algo- 
rithm in the Rhetenor module reconstructs the ancestral 
states along ES axes, plots all terminal and internal 
phylogenetic nodes into the morphospace, and connects 
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adjacent nodes by drawing branches between them. Phy- 
logenetic relationships were based on those reported by 
Rutschmann et al. (2011) and Matschiner et al. (2011). 
Branch lengths were calculated using mean value diver- 
gence dates reported by Matschiner et al. (2011). 

To assess whether disparity increases rapidly at an early 
stage in the icefish radiation and then asymptotes, as 
would be predicted in a scenario of rapid early diversifi- 
cation ("early burst") under conditions of ecological 
opportunity (Gavrilets and Losos 2009), we used DTT 
analyses to evaluate how shape disparity changed through 
time in comparison to trait evolution under a BM model. 
Analyses were implemented in R using the package "gei- 
ger" (Harmon et al. 2008) and the same phylogenetic 
framework as used for the phylomorphospace visualiza- 
tions. This method calculates disparity using average pair- 
wise Euclidean distances between species as a measure of 
variance in multivariate space (e.g., Zelditch et al. 2004). 
As input we used mean ES scores per species along axes 
ESI to ES8, encapsulating 95.8% of shape variance. Fol- 
lowing Harmon et al. (2003), relative disparities were cal- 
culated by dividing a subclade's disparity by the disparity 
of the entire clade. Relative subclade disparities were cal- 
culated for each node in the phylogeny, progressing up 
the tree from the root. At each node, the relative disparity 
value was calculated as the average of the relative dispari- 
ties of all subclades whose ancestral lineages were present 
at that time (Harmon et al. 2003: 961). Relative disparity 
values that are close to 0.0 indicate that subclades contain 
only a small proportion of the total variation and there- 
fore overlap in morphospace occupation is minimal 
between the different subclades, whereas, conversely, rela- 
tive disparity values that are close to 1.0 indicate extensive 
morphological overlap. To quantify how mean disparity 
compared to evolution under a BM model, 1000 simula- 
tions of morphological diversification were calculated on 
the phylogeny, and these theoretical subclade disparity 
values were plotted alongside the observed disparity 
values for opercle shape data. A morphological disparity 
index (MDI) metric was obtained, representing the area 
contained between the line connecting observed relative 
subclade disparity points versus the line connecting med- 
ian relative disparity points derived from BM simulations 
(Harmon et al. 2003). If the observed subclade disparity 
line plots above the BM line then the clades defined by 
that time slice have tended to generate higher disparity in 
the modern fauna than expected under the null and over- 
lap morphospace occupied by the overall clade. 

Model fitting 

BM, early burst (EB), and Ornstein-Uhlenbeck (OU) evo- 
lutionary models were fit to the data set of mean ESI 



scores for opercle shape. These models describe different 
processes of morphological evolution on a chosen phylog- 
eny and offer predictions about measures (e.g., disparity) 
of morphological trait evolution. The EB model predicts 
rapid morphological diversity early in the history of a 
group, followed by limited diversification as ecological 
niches are filled over time (e.g., Harmon et al. 2010). 
Under a BM model, trait evolution is simulated as a ran- 
dom walk and after each speciation event, the random 
walk continues independently of previous changes, and 
these changes are drawn from a normal distribution of 
zero and a variance proportional to branch length, hence 
phenotypic trait variance is predicted to increase with 
time in an unbounded fashion. The OU model is used to 
model stabilizing selection for a phenotypic trait value, 
and is similar to a BM model except traits are being 
pulled toward an optimal value, measured by a parameter 
(alpha) (Butler and King 2004; Hansen et al. 2008). 

Methods for modeling evolutionary processes are 
largely implementable only for univariate data and there- 
fore we chose ESI as representative of opercle shape 
because it represents the maximum variance in the sample 
(39.9%). We repeated model fitting also for ES2 (20.6%) 
to assess the consistency of the best chosen model. Akaike 
information criterion (AIC) values were used to compare 
the fit of each model to the data (Akaike 1974; Wagen- 
makers and Farrel 2004), and specifically we report a 
modified version, AICc, which performs better when the 
number of observations per parameter is small (Burnham 
and Anderson 2010; Hunt and Carrano 2010). The AICc 
values for each model were transformed into differences 
from the minimum observed AICc value A; 
(AICc) = AICc— min AICc. The differences were then 
transformed into AICc weights using the calculation: 



Wi(AICc) 



exp[-| x A ; (AICc)] 
E;exp[-|x A ; (AICc) 



The resulting values sum to one across a set of candi- 
date models, and can be interpreted as the proportional 
support received by each model (Hunt and Carrano 
2010). Model fitting was conducted using the function 
fitContinuousQ in the "geiger" package for R. 

Measurement error 

Error associated with the shape variables derived from 
outline data sets was calculated following the methodol- 
ogy of Arnqvist and Martensson (1998). Landmark data 
collection was replicated five times each for a subset of 
four specimens (A. nudiceps, A. skottsbergi, Chaenocepha- 
lus aceratus, and Dissostichus mawsoni), these were 
selected to include representatives from each of the four 
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families, and outlines were interpolated for the error 
repeats and added to the original data set. ES analysis was 
used to obtain shape variables and a one-way ANOVA 
was then performed on the outputted shape variables to 
detect whether the among-individual variance was greater 
than the within-individual (repeated) variance. The 
repeatability (R) value scales between 0 and 1. An R value 
of 0 would represent a sample in which all variance is 
found within individuals, whereas an R value of 1 would 
indicate all the variance is due to differences between 
individuals (see Wilson et al. 2011). 

Results 

Measurement error 

Measurement error was calculated across the first six ES 
axes (ES1-ES6) accounting for 91.8% of the total sample 
variance, and each comprising between 3% and 39.9% of 
variance. One-way ANOVAs conducted on a subsampled 
data set including all error replicates (N = 20) plus origi- 
nal outlines resulted in R values of between 0.90 and 
0.99, indicating a high level of replication for outline 
capture (Table S2). 

Patterns of opercle shape change 

The first three ES axes accounted for 75.3% of shape vari- 
ance in the sample. Shape variance along ESI (39.9%) 
was localized along two axes of the opercle outline. Nega- 
tive ESI scores reflected extension along a diagonal axis 
from the anterior dorsal margin to the posterior ventral 
margin of the bone coupled with compression along an 
axis from the posterior dorsal margin to the ventral tip. 
Conversely, positive ESI scores reflected compression 
along the anterior dorsal margin and posterior ventral 
margin, in addition to extension along the posterior 
dorsal margin and ventral tip (Fig. 3A). These differences 
resulted in separation between species belonging to Noto- 
theniidae, typically having negative scores along ESI, from 
members of Channichthyidae and Bathydraconidae, 
mostly characterized by positive ESI scores (Fig. 3A). 
Specifically, specimens of Notothenia rossii (Fig. 3A, label 
a) had the most extreme negative scores and specimens of 
C. aceratus the greatest positive scores along the axis 
(Fig. 3 A, label b). As for ESI, mean shape models for 
shape change along ES2, which represented 20.6% of 
shape variance in the sample, also indicated two alternat- 
ing axes of extension and compression along the opercle 
margin. Negative ES2 scores described extension along 
the entire dorsal margin of the opercle and lower portion 
of the ventral margin, alongside compression occurring 
broadly along the proximal margin and the upper portion 





Nototheniidae Artedidraconidae 
Bathydraconidae • Channichthyidae 



Figure 3. Phylomorphospace projections of notothenioid relationships 
on eigenshape (ES) axes ES1 and ES2 (A), and ES2 and ES3 (B) axes, 
describing interspecific differences in opercle shape. Branch lengths 
are taken from Matschiner et al. (2011), branches are colored by 
clade, and the root is denoted by concentric circles shaded black. 
Mean shape models illustrate, using vector displacements, the 
patterns of outline shape change associated with each axis. Tip labels, 
see Results for detail: a, Notothenia rossii; b, Chaenocephalus 
aceratus; c, Neopagetopsis ionah; d, Trematomus tokarevi; e, 
Trematomus eulepidotus. 

of the distal margin. Positive ES2 scores reflected changes 
along these axes in the opposite direction (i.e., compres- 
sion instead of extension, and vice versa). Similar to ESI, 
N. rossii also occupied the most negative portion of ES2, 
whereas specimens of Neopagetopsis ionah (Fig. 3A, label 
c) had the greatest positive scores, equating to a lateral 
extension of the distal tip of the operculum, resulting in a 
right-angled triangle shape appearance of the bone. ES3 
accounted for 14.9% of shape variance, and shape 
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differences included a combination of variance explained 
by ESI and ES2, thus resulting in two antagonistic modes 
of shape change occurring along each margin of the bone 
(Fig. 3B). 

Results from ANOVA tests performed on ES1-ES8 
scores, representing 95.8% of the sample variance, using 
"families" as groups indicated significant differences 
between Channichthyidae and Nototheniidae along ESI 
(^3,89 = 8.525, P < 0.001, Bonferroni corrected), ES2 
(^3,89 = 12.387, P < 0.001, Bonferroni corrected), and 
ES3 (_F 3>89 = 4.706, P < 0.001, Bonferroni corrected). 
Canonical variates analysis (CVA) performed on ES1-ES8 
scores using all specimens in the sample, resulted in three 
canonical functions that explained 100% of the sample 
variance. Only the first canonical function (eigen- 
value = 2.73) accounting for 95.6% of the variance was 
significant using Wilks' Lambda i'/ 2 ^ 89 = 119.46, 
P < 0.001) (Table S3). 

Disparity through time 

Phylomorphospace plots of ESI versus ES2 (Fig. 3A) and 
of ES2 versus ES3 (Fig. 3B) indicate a phylogenetic struc- 
turing of taxon distribution in shape space, particularly the 
separation of Nototheniidae and Channichthyidae and the 
distribution of Bathydraconidae and Artedidraconidae 
typically in-between those other two families. Average 
clade disparities for each clade were calculated from tip 
disparity values using the tip disparity function in the 
geiger package (per Harmon et al. 2003, 2008). These val- 
ues were summed for each of the four clades and shape 
disparity was found to be highest for the Nototheniidae 
(0.96), followed by the Channichthyidae (0.67), the Arte- 
didraconidae (0.16), and lastly the Bathydraconidae (0.11). 
Because sampling of species was unequal across the 
families, in part due to underlying differences in species 
diversity, the disparity values were subject to a simple stan- 
dardization by number of taxa in each clade to yield an 
average per species, which was highest for Channichthyidae 
(0.096), followed by Artedidraconidae (0.081), Notothenii- 
dae (0.074), and, lastly, Bathydraconiidae (0.055). 

The DTT method was used to assess how opercle shape 
and size disparity compared with expected disparity based 
on simulations using a neutral evolution BM model 
(Fig. 4). Overall, shape disparity using ES scores reflecting 
the positioning of taxa in multivariate shape space is 
greater than expected by BM simulations. A similar result 
is obtained using only size disparity. MDI values, calcu- 
lated as the area contained between the solid and dotted 
lines in Figure 4 or in other words the observed relative 
disparity points versus the line connecting median relative 
disparity points from the BM simulations, were similar 
for shape (0.341) and size data (0.453). 
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Figure 4. Disparity-through-time plot for opercle shape (solid black 
line) data, and opercle size from centroid size (solid red line) data. 
Mean values were used for species with more than one representative 
specimen. Disparity along the Y axis is the average subclade disparity 
divided by total clade disparity calculated at each internal node. The 
dotted line represents evolution of the data under Brownian motion 
(BM) simulations on the same phylogeny. Time values are relative time 
as per Harmon et al. (2003), whereby 0.0 represents the root and 1.0 
represents the tip. The most recent 20% of the plot was omitted to 
avoid the effect of "tip overdispersion" due to missing terminal taxa 
(Muschick et al. 2012). 



Evolutionary models 

The fit of the EB, OU, and BM models was assessed using 
the Akaike information criterion corrected for small sam- 
ple size (AICc), which can be used to compare models 
that have different numbers of parameters (BM has two 
parameters, OU has three) and therefore have noncompa- 
rable log likelihoods. AICc values indicate that the best fit 
to ESI shape data was the OU model (AICc = —23.02) 
followed by the BM model (AICc = -19.21) and lastly 
the EB model (AICc = -16.59) (Table 2). A similar result 
was found for ES2, also best supported by OU 
(AICc = -33.70), followed by BM (AICc = -21.69), and 
least supported by the EB model (AICc = —19.06). 
Results of AICc weight calculations indicated a compara- 
tively high probability that the OU model (0.84) was the 
best model given the data and the set of candidate models 
(Table 2). 

Patterns of shape change in relation to 
habitat and trophic niche inferred from 
stable isotope data 

A significant relationship was not found for results of 
PGLS regression analyses using stable isotope values for 
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Table 2. Comparison of evolutionary models fit to opercle shape 
data (ES1). Akaike weight was calculated from AlCc. 











Akaike 


Model 


AIC 


AlCc 


Log L 


weight 


Early Burst (EB) 


-17.79 


-16.59 


11.89 


0.034 


Brownian Motion (BM) 


-19.79 


-19.21 


11.90 


0.125 


Ornstein -Uhlenbeck (OU) 


-24.23 


-23.02 


15.11 


0.841 



<5 13 C and <5 15 N against the matrix of mean scores along 
ES1-ES8 for all species (r 2 < 0.15, P < 0.60). Members of 
the Channichthyidae and the Nototheniidae showed the 
greatest amount of spread along ESI and along <5 15 N 
values (Fig. 5 A) and a general, although not significant 
(P = 0.1493), trend of lower ESI scores associated with 
higher <5 15 N could be observed, indicating that species 
inferred to occupy higher trophic levels typically had 
opercles with elongated posterior portions of the dorsal 
margin and that tapered more sharply along the entire 
posterior margin (see Fig. 3 top-right mean shape model), 
although this was not evident for ES2 scores (Fig. 5B). 
Rutschmann et al. (2011) previously noted that species 
with lower d 13 C values were typically classified as pelagic, 
whereas benthic species were found to have higher <5 13 C 
values. Specific regions of morphospace were not exclu- 
sively occupied by benthic or pelagic species (Fig. 6). For 
instance, bathydraconids and artedidraconids are consid- 
ered the most benthic families within Notothenioidei 



(La Mesa et al. 2004), but occupied broadly average 
scores on ESI (Fig. 6A) and slightly higher than average 
scores on ES2 (Fig. 6B), although species with the highest 
ES2 scores occupied either a pelagic (N. ionah, Fig. 6B, 
label a) or benthopelagic niche (Cryodraco antarcticus, 
Fig. 6B, label b). Of note, C. aceratus, an exception 
among the largely pelagic Channichthyidae, is considered 
a benthic predator, mainly feeding on Champsocephalus 
gunnari (Reid et al. 2007), and is found to occupy sepa- 
rate regions of ESI (high positive score, Fig. 6A, label c) 
and ES2 (high negative score, Fig. 6B, label d) reflecting a 
slightly different opercle morphology to other members 
of the group. Labeling of specimens according to their 
feeding strategy indicates a broad overlap in opercle mor- 
phology between benthic and pelagic species, occupying 
mostly the area of -0.20 to 0.20 along ESI by -0.10 to 
0.10 along ES2 (Fig. 7). Semipelagic species, represented 
by Lepidonotothen larseni and N. rossii have low ESI and 
ES2 scores, forming a group slightly distinct from the 
benthic and pelagic species (Fig. 7) and equating to an 
opercle with an anterior margin tapering along its length 
in a posterior direction such that its most ventral tip is 
somewhat shifted posteriorly, compared to species with 
higher ES scores on these two axes. 

Discussion 

We investigated the evolution of opercle shape in the 
adaptive radiation of notothenioids by quantifying shape 
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Figure 5. Mean shape scores for each notothenioid species along eigenshape (ES) axes ES1 (A) and E2 (B) plotted against mean 5 15 N values, 
denoted per mil (% 0 ), taken from Rutschmann et al. (2011). Tip labels, see Results section for further detail: a, Neopagetopsis ionah; b, 
Chaenocephalus aceratus. 
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Figure 6. Mean shape scores for each notothenioid species along eigenshape (ES) axes ES1 (A) and E2 (B) plotted against mean <) 13 C values, 
denoted per mil (%„), taken from Rutschmann et al. (201 1). Tip labels, see Results section for further detail: a, Neopagetopsis ionah; b, Cryodraco 
antarcticus, c, d, Chaenocephalus aceratus. 



disparity, phylogenetic patterns of shape evolution, and 
ecological correlates in the form of stable isotope values 
to assess how ecological and morphological (shape) dis- 
parity are interrelated. Our focus on the evolutionary 
morphology of a craniofacial bone addresses how shape 
disparity data may inform our growing understanding 
of the features that define the adaptive radiation model 
or patterns that may be uncovered across different 
groups. 

Our main findings are that (1) DTT results show oper- 
cle shape and size disparity for subclades tended to gener- 
ate higher disparity in the modern fauna than would be 
expected under the neutral evolution BM model (Fig. 5), 
and evolutionary model comparisons indicate that the 
OU model is the best fit to our data and the "early burst" 
model is the least well supported, (2) the main evolution- 
ary axis of opercle shape change (ESI) reflects movement 
from a broad and rather more symmetrically tapered 
opercle to one that narrows along the distal margin, but 
with only a slight shape change on the proximal margin, 
(3) the distribution of taxa in shape space ordinations 
reveals a broad diversity of realizable opercle morphologi- 
es (Fig. 3) and phylomorphospace projections show clear 
phylogenetic groupings for opercle outline shape and a 
wide distribution of morphospace occupation for mem- 
bers of the family Nototheniidae, particularly extended by 
species belonging to the genus Notothenia, which occupy 
a portion of morphospace unexplored by other species 
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Figure 7. Plot of eigenshape (ES) axes ES1 and ES2 representing 
60.5% of the sample variance. Markers indicate feeding strategy 
taken from literature sources (Gon and Heemstra 1990; Reid et al. 
2007; Rutschmann et al. 2011). 



(Fig. 4), and (4) a significant relationship was not 
detected between opercle shape and isotope values using 
PGLS regression. 
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Opercle shape and benthic/pelagic trends 

In contrast to other morphological features that have 
been quantified in the classical examples of adaptive radi- 
ation such as cichlids and Anolis lizards, the study of evo- 
lutionary patterns of craniofacial bone shape has received 
comparatively less attention as previous studies have first 
focused on traits that are the likely candidates to display 
ecologically or functionally related variability, such as 
whole-body shape (Barluenga et al. 2006; Clabaut et al. 
2007; Berner et al. 2010; Harrod et al. 2010) or the jaw 
apparatus (Muschick et al. 2011, 2012). A notable excep- 
tion are the studies of Kimmel and others that have 
examined opercle variability (Kimmel et al. 2008; Arif 
et al. 2009; Kimmel et al. 2010) in different populations 
of three-spined sticklebacks (but see also Willacker et al. 

2010) , a well-established subject of study for speciation 
research (e.g., Schluter and McPhail 1992; Shapiro et al. 
2004; Colosimo et al. 2005). The major axis of shape vari- 
ation found in the opercle of three-spined stickleback 
populations from Iceland to diverse locations along the 
western coast of North America reflects a dilution- 
diminution mode of shape change (Kimmel et al. 2008, 

2011) , that is, an anterior-posterior stretching coupled 
with a dorsal-ventral compression of the outline shape. 
This pattern explains change between freshwater and mar- 
ine populations, whereas the second axis of shape change 
(PC2: Kimmel et al. 2011) is attributed to foraging ecol- 
ogy along the benthic-limnetic axis and translates to an 
overall widening of the opercle. Our mean shape models 
indicate that for notothenioids the major axis of shape 
variability (=ES1) in the sample reflects a similar exten- 
sion and compression, but these axes of shape change are 
not strictly in the craniocaudal and anterior-posterior 
direction, instead being slightly offset (Fig. 3). The gen- 
eral trend along ES2 also reflects a widening and narrow- 
ing of the opercle margin, as for sticklebacks (Kimmel 
et al. 2011). A lack of clear phylogenetic segregation in 
Figure 5A also indicates that along ESI members of the 
Channichthyidae and Nototheniidae therefore have 
evolved broadly similar opercle shapes in relation to their 
position along the pelagic-benthic axis (Fig. 6A). Besides 
sticklebacks, differences in feeding mechanism are already 
known to be reflected in body shape and bone morphol- 
ogy among benthic and limnetic morphotypes in cichlids 
(e.g., Barluenga et al. 2006; Clabaut et al. 2007; Muschick 
et al. 2012). The finding that benthic species in this study 
generally have an extended posterior margin of the oper- 
cle compared to pelagic species is consistent with the 
results of Klingenberg and Ekau (1996) who examined a 
series of body measurements among several Notothenii- 
dae belonging to the subfamilies Trematominae and Pleu- 
ragramminae. Klingenberg and Ekau (1996) found that 



benthic species had larger values for head width, which 
we here may consider to be reflected in the opercle by an 
extension of the posterior margin, and mouth length 
measures than pelagic species. Those authors speculated 
that these morphological features may reflect the larger 
sized prey available for consumption in benthic environ- 
ments. 

Evolutionary model fitting 

Our data indicate a strong preference for the OU model, 
which models selection to a single (global) optimum for 
all species, and suggests that the here observed disparity 
patterns may result from an adaptive peak or constraint, 
as highlighted more broadly in several other fish radia- 
tions, such as cichlids (Young et al. 2009; Cooper et al. 

2010) and in agreement with a recent broad-scale geomet- 
ric morphometric study of cranial and postcranial bone 
shape in actinopterygians (Sallan and Friedman 2012). 
Assuming that a single global optimum morphology is 
indeed accurate for notothenioids and given the benthic/ 
limnetic habitat variation in the clade (Rutschmann et al. 

2011) , one would not expect an association of opercle 
shape with habitat or diet, which is supported here by a 
lack of significant relationship between isotope values and 
opercle shape data. The OU model expects more evolu- 
tion to be apparent on later branches of phylogeny as 
selection to the optimum would result in phylogenetic 
signal generated from evolution at earlier branches being 
erased. Although the OU model supports the presence of 
an optimum, this conclusion must be taken cautiously 
here because the DTT results indicate disparity is concen- 
trated within subclades, that is, to say closely related 
species differ considerably in morphology. This conflicts 
with convergence to a single optimum (alpha), and hence 
we suggest support for the OU model may rather indicate 
loss of phylogenetic signal due to potentially rapid diver- 
gence rather than convergence to an optimum. 

At early stages of an adaptive radiation it is predicted 
under the "early burst" model that measures of disparity 
are high, followed by a subsequent drop in those values 
as time passes and available niche space falls to zero (e.g., 
Seehausen 2006; McPeek 2008). Model comparison results 
indicate that our data fit least well to this "early burst" 
model, which had the highest AICc value of all three 
models tested. Also, although we do find early peaks in 
opercle shape and size disparity (Fig. 4), which would be 
indicative of the rapid, early filling of empty niches, our 
plot does not support an "early burst" scenario (e.g., 
Gavrilets and Vose 2005) because we find a second peak 
in disparity occurring later in relative time (before 0.8, 
Fig. 4), and under an "early burst" scenario there 
would be little opportunity for subsequent ecological 
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diversification in subclades (Harmon et al. 2003; Burbrink 
and Pyron 2010). 

The second peak in disparity corresponds to the subc- 
lade within the family Nototheniidae including species of 
Trematomus, and the subclade comprising all representa- 
tive species of the Channichthyidae with the exception of 
Champsocephalus gunnari (Fig. 1). When examining the 
phylomorphospace plots for ESI and ES2 (Fig. 3A), mor- 
phospace occupation for the Channichthyidae is consider- 
ably extended by two taxa: N. ionah that displays low ESI 
values and high ES2 values (Fig. 3A, label c) and C. acera- 
tus that displays high ESI values and low ES2 values (top 
right of Fig. 3 A, label b). These two species may thus be 
contributing considerably to high values of disparity later 
in the DTT plot. Along with species of Notothenia, N. 
ionah also appears as an outlier on plots of i5 15 N versus 
ESI (Fig. 5A, label a), falling well below the majority of 
taxa in that plot. Similarly, the high score along ESI for C. 
aceratus, which as a top benthic predator (Kock 2005; Reid 
et al. 2007) stands out among the other largely pelagic 
channichthyids, results in that species being located out- 
side (above) the main group in Figure 5A (label b). In the 
case of Trematomus, here represented by six species, 
Rutschmann et al. (2011) showed that species of this 
genus were differentiated in isotopic signatures, indicating 
trophic niche separation within the genus or a large niche 
space, and reports of stomach contents for different spe- 
cies corroborate this finding (Brenner et al. 2001). Within 
our sample, the phylomorphospace plot indicates consid- 
erable variation particularly in ES2 scores among members 
of Trematomus, especially T. tokarevi (benthic, Fig. 3A 
label d) compared to T. eulepidotus (epibenthic/pelagic, 
Fig. 3 A label e), and these differences may have contrib- 
uted to elevated disparity for that node. Near et al. (2012) 
conducted a series of DTT analyses on buoyancy measures 
for 54 species of notothenioids and similarly their plots 
(Near et al. 2012: Fig. 3A-C) also revealed a second peak 
in disparity, particularly for Channichthyidae and species 
of Trematomus, which those authors related to the 
repeated colonization of benthic, epibenthic, semipelagic, 
and pelagic habitats among closely related lineages. The 
latter is thought to have happened as a consequence of the 
repeated creation of open niches following extinctions 
caused by icebergs and glaciers scouring the continental 
shelf and decimating near-shore fauna (Tripati et al. 2009; 
Near et al. 2012). 

More broadly, the lack of an "early burst" pattern in 
our data set fits with the results of Harmon et al. (2010), 
who performed a broad survey of 49 animal clades and 
found little evidence of an "early burst" model of mor- 
phological change, and recently Ingram et al. (2012) sug- 
gested that this may be explained by the ubiquity of 
omnivory in natural food webs. Ingram et al. (2012) 



found that the "early burst" scenario was not detected for 
clades containing many omnivorous species that fed at 
multiple trophic levels; a feature common also for noto- 
thenioids, which include several species that feed oppor- 
tunistically throughout the water column (e.g., Eastman 
2005). Although omnivory was suggested as one possible 
determinant of the adaptive burst scenario, a general 
trend hinted by those results is that the persistence of an 
"early burst" pattern may be related to the relative extent 
to which niche axes (such as diet, microhabitat, and cli- 
mate) are distinct and stable over time (Ingram et al. 
2012). 

Patterns of diversification in notothenioids 

The constituent groups of the notothenioid radiation 
have undergone different amounts of ecological and mor- 
phological diversification, with some, such as the artedi- 
draconids that are all sedentary benthic fishes, displaying 
little (Eastman 2005). Our disparity values and phylomor- 
phospace plots to some extent reflect these patterns, 
particularly for the notothenioids, which display the high- 
est disparity values and the most expanded occupation of 
morphospace (Fig. 3). Notothenioids are ecologically 
diverse and include benthic (around 50% of within-group 
species diversity, Eastman 1993), epibenthic, semipelagic, 
cryopelagic, and pelagic forms. They are also the only 
group containing species that have so far been determined 
as neutrally buoyant (Pleuragramma antarcticum and 
D. mawsoni are examples in our study), a feature that has 
been achieved, despite not possessing a swim bladder, 
through reduced skeletal mineralization and lipid deposi- 
tion (DeVries and Eastman 1978; Eastman and DeVries 
1982; Eastman 1993). Most distinct in our morphospace 
plots is the location of Notothenia species that typically 
have an opercle that widens at the posterior margin (ESI) 
and has a posteroventrally tapering dorsal margin (see 
top-left mean shape model, Fig. 3A). Representing the 
opposite end of the body mass scale compared to the 
neutrally buoyant members of the Nototheniidae, species 
of Notothenia are large, heavy fishes that are able to move 
up and down in the water column to feed on both pelagic 
and benthic prey, and are able to alter their diet in rela- 
tion to prey availability (e.g., Fanta et al. 2003). Notothe- 
nia coriiceps, for example, is known to feed on 
macroalgae, most likely to ingest also the associated 
amphipods more efficiently (Iken et al. 1997; Fanta et al. 
2003), when its preferred food source of krill is unavail- 
able. Notothenia rossii also ingests different food during 
its juvenile stages, switching from a pelagic to largely ben- 
thic habit in adulthood, which may have further implica- 
tions for opercle and craniofacial bone development in 
general. Burchett (1983) examined this ontogenetic shift 



2013 The Authors. Ecology and Evolution published by John Wiley & Sons Ltd. 



3177 



Opercular Shape Evolution in Icefishes 



L. A. B. Wilson ef al. 



from pelagic to benthic lifestyle and found an associated 
change in head shape (length and diameter) and a deep- 
ening of the body over the course of ontogeny. The main 
result of the foraging habit versus opercle shape plot, 
showing broad overlap in opercle morphology among 
different foraging categories (Fig. 7), is perhaps not 
unsurprising, given the dietary plasticity of many notothe- 
nioids (Eastman 2005), the aforementioned Notothenia 
being an excellent example (e.g., Foster and Montgomery 
1993). The most logical reasoning behind the range of 
morphotypes is that notothenioids inhabit an ecosystem 
with relatively low species diversity and reduced competi- 
tion, both of which would not act to accelerate ecomor- 
phological divergence (Eastman 2005) to the degree 
found among other radiations. 

Conclusions 

A major impetus for the study of adaptive radiations is to 
uncover generalized patterns among different groups. In 
this way, common features may speak for the importance 
of a given process in the generation of morphological 
diversity (Gavrilets and Losos 2009). Here, we use out- 
line-based geometric morphometries to quantify opercle 
shape across notothenioids. We identify axes of shape 
change, particularly a widening of the opercle bone, that 
have been recovered in other adaptive radiations (three- 
spined sticklebacks) and a trend in opercle shape change 
along the benthic-pelagic axis, underlining the impor- 
tance of this axis for diversification in notothenioids. We 
find that opercle shape and size disparity for subclades 
tended to generate higher disparity in the modern fauna 
than would be expected under neutral evolution, and that 
the OU model best fits the evolution of opercle shape. 
Support for the OU model may reflect loss of phyloge- 
netic signal due to potentially rapid divergence. Opercle 
shape represents one of few features that can be quantita- 
tively assessed for both extant and extinct species flocks 
(Wilson et al. 2013b), and therefore provides an especially 
useful opportunity for integrative study between evolu- 
tionary biology and paleontology (e.g., Sanchez- Villagra 
2010; Wilson 2013b), an approach that has yet to be fully 
explored in the context of adaptive radiation, and one 
that holds potential to yield valuable insights into modes 
of species diversification in deep time. 
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